#ifndef _SYSTEM_M_H_
#define _SYSTEM_M_H_
#include<string>
#include"help.h"
#include"base.h"
#include"orbital.h"
#include"twobodystate.h"
#include"system_ho.h"

using std::string;

template<class OrbitalType,class DataType>
  class System_TzM:public System<OrbitalType,TwoBodyState_TzParM,DataType>
{
 public:
  typedef TwoBodyState_TzParM TwoBodyStateType;
  typedef TwoBodyMatElement<DataType> TwoBodyMatEle;
  
  void setupTwoBodyStates();

  /// setup IndiceIn2BStates.
  void setupIndiceIn2BStates();
  vector<int> IndiceIn2BStates;

  void setupTwoBodyMatSketch();
  
  
  void get2BmatAt(int a,int b,int c,int d,TwoBodyMatEle & MatEle) const;
  void set2BmatAt(int a,int b,int c,int d,const TwoBodyMatEle & MatEle);

 private:
  int mapab(int a,int b) const//a<b,a map from a<b matrix to one dimension vector
  {
    int totalOrbitals=this->Orbitals.size();
    return a*totalOrbitals+b-1-(a*(a+3))/2;
  }
};

template<class OrbitalType,class DataType>
  void System_TzM<OrbitalType,DataType>::setupTwoBodyStates()
{
  this->TwoBodyStates.clear();
  for(int i=0;i<this->Orbitals.size();i++)
    for(int j=i+1;j<this->Orbitals.size();j++)
      {
	int Tz=(this->Orbitals[i].tz + this->Orbitals[j].tz)/2;
	int Par=(this->Orbitals[i].par + this->Orbitals[j].par)%2;
	int M=(this->Orbitals[i].mm + this->Orbitals[j].mm)/2;;
	this->TwoBodyStates.push_back(TwoBodyStateType(i,j,Tz,Par,M));
      }
  this->setupTwoBodyChannels();//noting that the twobody states are sorted after this calling.
}

template<class OrbitalType,class DataType>
  void System_TzM<OrbitalType,DataType>::setupIndiceIn2BStates()
{
  int totalOrbitals=this->Orbitals.size();
  IndiceIn2BStates.resize((totalOrbitals*(totalOrbitals-1))/2);
  for(int i=0;i<this->TwoBodyStates.size();i++)
    {
      int a=this->TwoBodyStates[i].a;
      int b=this->TwoBodyStates[i].b;
      //a<b
      IndiceIn2BStates[mapab(a,b)]=i;
    }
}

template<class OrbitalType,class DataType>
void System_TzM<OrbitalType,DataType>::get2BmatAt(int a,int b,int c,int d,TwoBodyMatEle & MatEle) const
{
  MatEle.setZero();
  if(a==b || c==d) return;//pauli principle
  int phase=1;
  //make sure a<b,c<d,(ab)<=(cd)
  if(a>b)
    {
      swap(a,b);
      phase*=-1;
    }
  if(c>d)
    {
      swap(c,d);
      phase*=-1;
    }
  int Mab=(this->Orbitals[a].mm + this->Orbitals[b].mm)/2;
  int Mcd=(this->Orbitals[c].mm + this->Orbitals[d].mm)/2;

  if( Mab!=Mcd ) return;
  int ab=IndiceIn2BStates[mapab(a,b)];
  int cd=IndiceIn2BStates[mapab(c,d)];
  if(ab==-1||cd==-1) return;
  int channelab=this->ChannelIndice[ab];
  int channelcd=this->ChannelIndice[cd];
  //channelab should = channelcd
  if(channelab!=channelcd) return;
  int bra=this->TwoBodyStateIndice[ab];
  int ket=this->TwoBodyStateIndice[cd];
  if(bra>ket)
    swap(bra,ket);
  this->get2Bmat(channelab,bra,ket,MatEle);
  MatEle*=phase;
}

template<class OrbitalType,class DataType>
void System_TzM<OrbitalType,DataType>::set2BmatAt(int a,int b,int c,int d,const TwoBodyMatEle & MatEle)
{
  if(a==b || c==d) return;//pauli principle
  int phase=1;
  //make sure a<b,c<d,(ab)<=(cd)
  if(a>b)
    {
      swap(a,b);
      phase*=-1;
    }
  if(c>d)
    {
      swap(c,d);
      phase*=-1;
    }
  int Mab=(this->Orbitals[a].mm + this->Orbitals[b].mm)/2;
  int Mcd=(this->Orbitals[c].mm + this->Orbitals[d].mm)/2;

  if( Mab!=Mcd ) return;
  int ab=IndiceIn2BStates[mapab(a,b)];
  int cd=IndiceIn2BStates[mapab(c,d)];
  if(ab==-1||cd==-1) return;
  int channelab=this->ChannelIndice[ab];
  int channelcd=this->ChannelIndice[cd];
  //channelab should = channelcd
  if(channelab!=channelcd) return;
  int bra=this->TwoBodyStateIndice[ab];
  int ket=this->TwoBodyStateIndice[cd];
  if(bra>ket)
    swap(bra,ket);
  TwoBodyMatEle MatEleTemp=MatEle;
  MatEleTemp*=phase;  
  this->set2Bmat(channelab,bra,ket,MatEleTemp);
}


template<class OrbitalType,class DataType>
  void System_TzM<OrbitalType,DataType>::setupTwoBodyMatSketch()
{
  setupIndiceIn2BStates();
  int totalChannels=this->TwoBodyChannels.size();
  this->TwoBodyMat.resize(totalChannels);
  for(int channel=0;channel<totalChannels;channel++)
    {
      int num2BInChannel=this->TwoBodyChannels[channel].size();
      this->TwoBodyMat[channel].clear();
      this->TwoBodyMat[channel].resize( (num2BInChannel*(num2BInChannel+1))/2 );
    }
}



////****************************************************************
////****************************************************************

  ///orbitals with par,mm,tz
class Orbital_TzM:public Orbital_parmmtz
{
 public:
  typedef double DataType;  
  Orbital_TzM(){}
 Orbital_TzM(int _par,int _mm,int _tz,DataType _e):Orbital_parmmtz(_par,_mm,_tz),e(_e){}

  DataType e;
};

  
  
  ///system hold interaction in M scheme
class System_TzM_HO:public System_TzM<Orbital_TzM,double>
{
 public:
  typedef double DataType;
  typedef System_Tz_HO VSystem;//interaction system
  typedef Orbital_TzM OrbitalType;
  typedef TwoBodyMatElement<DataType> TwoBodyMatEle;
  const static int MatEledim;

 System_TzM_HO(VSystem *Ptr):pSystem(Ptr),hbar_omega(Ptr->hbar_omega){}
  System_TzM_HO(){}

  void setup(VSystem *Ptr)
  {
    pSystem=Ptr;
    hbar_omega=Ptr->hbar_omega;
    setup();
  }
  
  void setup()
  {
    setupOrbitals();
    setupTwoBodyStates();
    setupTwoBodyMatSketch();
    //setupTwoBodyMat();
  }
  void setupOrbitals();
  void setupTwoBodyMat();
    
  
  void get2BmatOnNeed(int a,int b,int c,int d,TwoBodyMatEle & MatEle) const;
  
  DataType get2BV(int A,int a,int b,int c,int d)
  {
    TwoBodyMatEle MatEle;
    get2BmatAt(a,b,c,d,MatEle);//try fetch directly
    if(MatEle.empty())//if it's empty, maybe it's not setup yet
      {
	get2BmatOnNeed(a,b,c,d,MatEle);
	set2BmatAt(a,b,c,d,MatEle);//setup
      }
    if(MatEle.empty()) return 0;
    return MatEle[0]-MatEle[3]*pSystem->hbar_omega/A;
  }

  DataType get1B(int A,int a,int b) const
  {
    if(Orbitals[a].isInAGroupWith(Orbitals[b]))
      return pSystem->get1B(A,Orders[a],Orders[b]);
    else
      return 0;
  }

  double get1B_p2(int bra,int ket) const;
  double get1B_p2_plus_r2(int bra,int ket) const
  {
    return bra==ket?Orbitals[bra].e:0;
  }


  void print2B() const;
  
  double jplus(int a,int b) const;///< J+ matrix element between two orbital a and b
  double jminus(int a,int b) const;///< J+ matrix element between two orbital a and b

  DataType rL(int bra,int ket,int L);
  
  void printOrbitals(int num=-1) const;
  void printTwoBodyStates() const;
      
  VSystem * pSystem;
  vector<int> Orders;//split from which orbital in pSystem->Orbitals, the main shell of each orbital

  double hbar_omega;
 private:
};

#endif
